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Abstract 



We present a multigrid algorithm for self consistent solution of the Kohn- 
Sham equations in real space. The entire problem is discretized on a real 
space mesh with a high order finite difference representation. The resulting 
self consistent equations are solved on a heirarchy of grids of increasing reso- 
lution with a nonlinear Full Approximation Scheme, Full Multigrid algorithm. 
The self consistency is effected by updates of the Poisson equation and the 
exchange correlation potential at the end of each eigenfunction correction cy- 
cle. The algorithm leads to highly efficient solution of the equations, whereby 
the ground state electron distribution is obtained in only two or three self 
consistency iterations on the finest scale. 

PACS numbers:31.15.Ar, 71.15.-m, 71.15.Nc, 02.70.-c 
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The plane wave pseudopotential method has proven to be an excellent computational 
strategy for solving large scale electronic structure problems in condensed phases No- 
table strengths of the method are the ability to use fast Fourier transforms for updating the 
self consistency equations, lack of dependence of the basis on atom positions, and the clear 
control of convergence with the cutoff energy determined by the shortest wavelength mode. 
However, the method can encounter difficulties in treating widely varying length scales. The 
problems are especially severe for surfaces, clusters, or systems where the pseudopotential 
varies rapidly near the nucleus such as first row elements or transition metals. Also, charged 
systems present significant complications in a plane wave code. In recent years, consider- 
able effort has been focused on alternative real space approaches utilizing finite elements or 
finite difference representations iPHlOU- Advantages of these approaches include the ability 
to handle finite or periodic systems with equal effort and the locality of each iteration step. 
Locality leads to simplicity in developing domain decomposition parallel algorithms. In ad- 
dition, it is relatively straightforward to implement adaptive grid refinement techniques in 
order to focus effort in spacial regions with large variations in the computed functions, for 
example near the nuclei in electronic structure computations [§],|11|. Finally, representation 
directly in real space allows for the use of multigrid (MG) algorithms with their excellent 
convergence characteristics and scaling properties 

Several groups have employed MG algorithms for various portions of the solution process 
for self consistent electronic structure computations. White, et al. developed a finite ele- 
ments method for electronic structure which utilized an MG solver for the Poisson equation. 
Davstad discretized the Hartree-Fock equations for diatomic molecules and solved the 
resulting two dimensional equations with a combination of the MG method and a Krylov 
subspace method for the coarsest grid equations. Gygi and Galli proposed an adaptive 
coordinate approach which places increased resolution near the nuclei using curved grids. 
They solved the Poisson portion of the problem with MG techniques. Briggs, et al. have 
developed an MG solver for the Kohn-Sham (KS) equations which utilizes MG methods 
for both the eigenvalue and Poisson problems. Their MG eigenvalue solver uses a double 
discretization approach and a linearized version of the MG method. Recently, Ancilotto, et 
al. have presented a similar method for solution of the KS equations, and have applied 
the method to calculations on charged lithium clusters. They show that the MG approach 
is as accurate and more efficient than the corresponding Car-Parrinello method. Modine, 
et al. 0] developed an adaptive grid real space method which employs MG preconditioning 
for solution of the Poisson equation. Finally, we have developed high order MG methods 
for solving the KS equations including all the electrons |10|. Typical of the MG solvers 



to date is the requirement of 20 or more self consistency cycles to obtain convergence. In 
this letter, we present a high order multigrid algorithm for solving the self consistent Kohn- 
Sham equations which obtains convergence in an order of magnitude less numerical effort. 
The method utilizes the nonlinear Full Approximation Scheme (FAS), full multigrid (FMG) 
eigenvalue technique of Brandt, et al. with the further inclusion of new methods for the 
self consistency portion of the problem. 

The first step in development of the MG algorithm is the discretization of the problem in 
real space. Here we employ a high order finite difference representation. This approach has 
been shown to yield accurate results in pseudopotential and all electron calculations P,p!0|. 
Consider the KS equations ITBI in the Local Density Approximation (LDA): 
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[-\v^ + v,jf{v)]A{v) = eM^). (1) 
The one-electron effective potential is: 

= Vesij) + V^^{p{j)), (2) 

with the electrostatic portion of the potential (nuclear potential plus Hartree potential) 
given by: 

I \ * I I I 
This potential can be obtained by real space numerical solution of the Poisson equation: 

V\es{r) = -Aivptotir), (4) 

where ptot is the total charge density on the grid due to the electrons and nuclei. The 
exchange- correlation potential in our spin restricted LDA computations is calculated with 



the VWN functional \U 



The Laplacian operator is discretized with a high order finite difference representation 
15| . In the present work, a 12th order form is used. The high order expression yields 
a large gain in accuracy, which reduces the number of required grid points significantly. 
The same Laplacian is used for solution of the Poisson and KS equations, resulting in a 
consistent level of accuracy throughout the solver. Due to this consistent representation 
on all levels, the only parameter controlling the accuracy of the solution is the fine grid 
spacing h. The other quantities are diagonal in the coordinate representation. The real 
space discretization scheme offers efficient calculation of electron densities and, as a result, 
the exchange-correlation potentials. The nuclear charge density is the discretized form of 
the Dirac delta function. The wavefunctions are vectors of length Ng, the total number 
of grid points on a given level. Integration is performed by simple trapezoidal summation 
on the three dimensional domains. The boundary conditions on the orbitals for the finite 
systems examined here are set to zero, while the electrostatic potentials on the boundary 
are obtained via multipole expansion of the total charge density to quadrupole order. 

The KS eigenvalue problem is nonlinear since one simultaneously solves for both the 
eigenvalues and the eigenfunctions. Therefore, we have employed the FAS eigenvalue tech- 
nique of Brandt, et al. This method allows for solution of nonlinear problems with 



efficiencies similar to linear ones. Here we outline the FAS method of solution for the eigen- 
value and Poisson problems. Both solvers are required for the self consistent problem. 
The discrete equations can be represented on the current finest grid as: 

L^U^ = f. (5) 

For the Poisson problem is the finite difference Laplacian on the fine grid with spacing 
h, is the electrostatic potential t'es(r), and f'^ is — 47rptot(r). Upper case for the potential 
denotes the exact numerical solution. Below, lower case implies the current approximation. 
In the eigenvalue equations becomes Uj^, the eigenfunctions, where the index i indicates 
the orbital number. The eigenvalue operator is the Hamiltonian minus Aj, where Aj is 
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the eigenvalue for orbital i, and there is no source term / . The eigenvalue has the same 
value on all levels at convergence and therefore is not labeled with a grid size. 

In the FAS method, the desired functions are represented on each grid level. The discrete 
equations on the coarse level, however, are modified by the inclusion of a defect correction 
term which is required to yield zero correction at convergence. Suppose we have an ap- 
proximate solution on the fine grid h, u^, which is obtained by relaxation sweeps on the 
grid: 

L'^u'' = /^ (6) 

We have used Gauss-Seidel or Successive Over- Relaxation updates for these smoothing steps. 
The coarse grid problem with grid spacing H = 2h can then be constructed as: 

LV^ll^f + r", (7) 

where J/f is the restriction operator which takes a local average of the fine grid function. Full 
weighting restriction is employed here, which weights the local points according to trapezoid 
rule integration. The defect correction is: 

= L^'l^u'^ - I^L^uK (8) 

The coarse grid equation is iterated to obtain an improved solution at that level. In a two 
grid method, the fine grid function is then corrected via: 

u^^u^-rI^H{u''-I^u% (9) 

where is the interpolation operator. These operations can be recursively extended to yet 
coarser grids. By inclusion of multiple levels, errors of all wavelengths are rapidly damped 
in the fine grid function. The only change required on grids two or further from the finest 
grid is that an additional term //f r'^ must be included in the defect correction to incorporate 
information from the previous level. 

In solving the eigenvalue equations in KS theory, some additional components are re- 
quired. On all grids except the coarsest, the orbital equations are updated just as outlined 
above. However, on the coarsest level, constraints must be imposed to maintain wavefunc- 
tion separation. If one had the exact solution on a fine grid and then restricted the 
orbitals to the next coarser level, the functions would no longer be orthogonal. Hence, the 
constraints on the coarsest level can be implemented by solving a linear matrix equation: 

{ufj^u^,)^{I^ull^u)) (10) 

These equations, since they are implemented only on the coarsest scale, require very little 
computational overhead, and can be solved either via exact matrix inversion or by iterative 
methods. The matrix size is of order g x g, where q is the number of orbitals. The eigenvalue 
is also updated on the coarsest level by computing the Rayleigh quotient: 



Finally, the convergence can be improved by including subspace orthogonalization via 
Ritz projection. At the end of an MG correction cycle, a Gram-Schmidt orthogonalization 
is first performed on the finest level, followed by diagonalization of the Hamiltonian in the 
basis of the occupied orbitals: 

u^LuoZi — XiZi = (12) 

where u symbolizes the q x Ng matrix of eigenf unctions on the grid, and the Zi are the 
coefficients used to improve the occupied subspace. This matrix problem is thus also of size 
q X q. 

The FMG approach is initiated by first iterating the self consistent problem on the coars- 
est level used for the eigenvalue problem. The initial approximation for the eigenf unctions 
is a set of random numbers of magnitude one. The coarsest level is chosen so that the 
eigenvalues have the same ordering as on the finest level. For the systems studied here, a 
coarse grid of 8^ or 17^ is required, depending on the problem. The Poisson solver extends 
the levels up to a 3^ grid where only the one central point is iterated. In most of these 
calculations the finest grid employed is of size 65'^, so the eigenvalue solver comprises three 
or four levels, while a total of six levels is utilized for the Poisson solver on the finest scale. 
During the initial iterations on the coarest level, the eigenfunctions are first relaxed and 
orthogonalized via a Gram-Schmidt process, and then the potential is updated. Once an 
initial approximation is obtained on the coarse level, the solution is interpolated to the next 
finer level, where the MG correction cycles are initiated. 

At the beginning of each MG V-cycle, the effective potential is updated once upon entry 
to the new finest level with the FAS method. Subsequently, the potential is updated at the 
end of each eigenfunction correction cycle; since the changes to the potential are relatively 
smooth, the potential corrections are performed with a simple V-cycle using the previous 
potential as input. A schematic diagram of the MG solver is presented in Figure 1. The 
Poisson equation is typically solved with only a few relaxation sweeps on the finest scale 
on either side of the V-cycle. During each eigenfunction correction cycle, the eigenfunctions 
are relaxed 3 times on each side of the V cycle. The solution of the entire electrostatic 
portion of the problem thus requires similar computational effort as for the update of a 
single eigenfunction. This algorithm differs from those of Refs. |^,^ by inclusion of updates 
of the eigenvalues and enforcement of the constraints on the coarsest level only without 
resorting to linearization of the equations. 

We have explored two approaches for the eigenfunction correction cycles. In the first, each 
eigenfunction is carried through the V cycle starting with the lowest energy function, and the 
effective potential is updated at the conclusion of the cycle. This sequential method leads 
to a rapidly convergent effective potential since the low lying states are quickly stabilized. 
In the second approach, all eigenfunctions are corrected simultaneously, and the effective 
potential is updated upon conclusion of the V cycle. For small systems the first approach 
is preferred since the Poisson updates are quite inexpensive. However, as we discuss below, 
the scaling of the two approaches differs, and the second method is preferable for larger 
systems. Its convergence behavior is slightly less dramatic, requiring three or four self 
consistency updates as opposed to two for the first method. The algorithm can readily be 
adapted for either approach depending on the problem. 

We have performed calculations on atoms, ions, and small molecules to test the conver- 
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gence and accuracy of the method. The numerical results presented here were obtained with 
the second approach discussed above. The convergence behavior is illustrated in Fig. 2 in 
calculations on the 4 electron Be atom and the 14 electron CO molecule. Comparison is made 
to recent MG and Car-Parrinello pseudopotential calculations |§] on the 8 valence electron 
C2 molecule. The FAS-FMG approach leads to significant acceleration; we emphasize that 
the present calculations include the nuclear singularity in the effective potential, making it a 
more challenging scenario for convergence to the ground state. We computed total energies 
for atoms and ions to obtain atomic ionization potentials at the LDA-VWN level (Table I). 
The total energy obtained for the He atom is -2.8345 au vs. the exact value of -2.8348 au. 
The computations on the ions are performed with equal effort as for the neutral species, and 
the IP results are of satisfactory accuracy [|l^. Finally, we present the eigenvalues of the 
CO molecule computed at the Xa level for comparison with previous numerical work |17 



(Table I). The computation was performed on a large 129^ domain since the dipole moment 
of CO is a sensitive function of the overall domain size. We obtained a dipole value of 0.25 
D C^0+ compared to the previous result of 0.24 D in fully converged numerical calculations 
p!8| , and 0.10 D in finite difference pseudopotential calculations [Q. As discussed in Ref. [Q, 
the most severe errors in real space all electron calculations occur in the regions around the 
nuclei. The limitations on the accuracy of our uniform domain results is a consequence of 
these errors. Two improvements are suggested: inclusion of pseudopotential techniques to 
remove the core and/or local grid refinements in the neighborhood of the nucleus [§,|ll|]- 
have generalized a second order multigrid local mesh refinement method [|1^] to arbitrary 



order and will include these refinements in the KS solver in future work ||TT|. Since the 
refinements are truly local, they require only modest computational overhead. 

To conclude, we discuss the scaling properties of the various steps in the FAS-FMG 
algorithm for the case where the orbitals span the whole domain. With q orbitals and A^^ 
fine grid points in the domain, the scalings of the important components of the algorithm are 
as follows: relaxation of orbitals {qNg), relaxation of potential (Ng), Gram-Schmidt process 
(q'^Ng), Ritz projection {q^Ng to construct the matrix and q^ to solve), computation of the 
eigenvalues on the coarse grid {qN^), and solution of the constraint equations on the coarse 
grid {q^N^ to construct the constraint matrix and q^ to solve). The most costly portion of 
the algorithm for the small systems examined here is the relaxation step for the orbitals. 
For large systems where orbital localization is possible, each of the steps becomes linear 
scaling with the number of electrons. The two algorithms for implementing self consistency 
discussed above differ in that the first approach, while exhibiting stronger convergence to 
the ground state, leads to a qNg scaling which persists even with localized orbitals. 

We have presented a new method for solution of the KS equations in real space which 
utilizes nonlinear multigrid techniques to solve the self consistent equations. The purpose 
has been to illustrate the rapid convergence behavior of the FAS-FMG algorithm in relation 
to other electronic structure methods. The efficiency is a result of the preconditioning on 
coarser levels and nonlinear treatment during the MG correction cycles. The combination 
of real space discretization with FAS-FMG multiscale techniques of solution leads to order 
of magnitude improvement in convergence behavior relative to other MG methods and the 
Car-Parrinello approach, and it should thus prove a useful numerical stategy for solving 
large scale electronic structure problems. 

We would like to thank Prof. Achi Brandt for many helpful discussions concerning this 
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TABLES 





h Ref. 16 


FAS-FMG 


He^ He+ 


0.15 0.969 


0.974 


Li^ Li+ 


0.35 0.198 


0.193 


0^ 02+ 


0.12 1.843 


1.821 




Orbital 


Ref. 17 


FAS-FMG 


ls(0) 


-18.745 


-18.832 


ls(C) 


-9.912 


-9.940 


a{2s) 


-1.045 


-1.060 


a* {2s) 


-0.489 


-0.494 


7r(2p) 


-0.413 


-0.407 


(7(2p) 


-0.304 


-0.302 


TABLE I. Calculations on Atomic Ionization Potentials and Eigenvalues for the CO molecule. 


The LDA calculations of Ref. 16 used a form for the correlation energy which interpolated between 


the Wigner 


form for low densities and the Gell-Mann/Brueckner form at high densities. For the 


CO calculation, the grid spacing was h = 0.1S35au and the bond length was taken 


as 1.13A. All 



energies are in hartrees. 
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FIGURES 
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FIG. 1. The FMG scheme. The fine grid solution is initiated with the interpolated functions 
from the coarse grid approximation(bottom of figure). At the end of each V cycle the potential is 
updated (indicated by x). 
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FIG. 2. Convergence behavior. The top curve is the Car-Parrinello result of Ref. 8. The second 
curve is the MG result of Ref. 8. The next is our result for the CO molecule with the FAS-FMG 
solver. The bottom curve is the FAS-FMG result for the Be atom. 
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